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Abstract 

A comprehensive approach for real-time computations using a database of parameterized linear reduced- 
order models (ROMs) is proposed. The method proceeds by sampling offline ROMs for specific values 
of the parameters and interpolating online the associated reduced operators. In the offline phase, a pre¬ 
processing step transforms the reduced operators into consistent sets of generalized coordinates prior to their 
interpolation. The present paper also introduces a consistency enforcement approach for models defined on 
arbitrary underlying meshes. In the online phase, the operators are interpolated on matrix manifolds. The 
proposed framework is illustrated on two realistic multi-physics problems: an inverse acoustic scattering 
problem around a submarine and flutter predictions for a wing-tank system. The second application is 
implemented on a mobile device, illustrating the capability of the proposed framework to operate in real¬ 
time. 

Keywords: Parametric model order reduction, database, interpolation, mobile computing, aeroelasticity, 
acoustic scattering 


1. Introduction 

Many engineering applications require the ability to generate predictions of the behavior of physical 
systems in real-time. Among those applications, one can mention design optimization, optimal control, 
the solution of inverse problems as well as uncertainty quantification. All of these applications require a 
large number of predictions for varying values of operating conditions. The operating conditions, usually 
described by a set of parameters, may define boundary conditions, initial conditions, physical or shape 
parameters that define the problem of interest and its underlying differential equations. However, each of 
these predictions usually demands computationally intensive calculations as accurate discretization of the 
underlying differential equations often leads to large scale systems of equations. 

Projection-based model reduction mm reduces the large computational cost associated with each solution 
of the underlying high-dimensional model (HDM) by reducing the number of degrees of freedom in the 
computation. For this purpose, a reduced-order basis (ROB) is defined and the solution is restricted to 
the subspace described by the ROB. The current most challenging model reduction problems are those 
associated with nonlinear systems and parameter variations. Nonlinear systems require additional levels of 
approximation to enable large computational speedups El a 0161 El m. The model reduction of parameterized 
systems is also challenging due to the non-robustness of reduced-order model with respect to parameter 
variations that requires an appropriate offline training phase 0 nni m Eg. Approaches addressing the 
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model reduction of nonlinear parameterized systems are proposed in [m na na HE]. The focus of the 
present paper is on the efficient model reduction of linear parameterized systems. 

More specifically, for linear parameterized systems, database approaches can be developed by pre- 
computing in an offline phase the reduced linear operators of a common reduced dimension for specific 
values of the parameters [ini Ha HU HD Ha uni ED •These linear operators are subsequently interpolated in 
the online phase for values of the parameters not present in the database. The small dimensionality of the 
reduced operators leads to real-time interpolation and predictions on-the-fly. 

The interpolation of local reduced operators is however a challenging because each reduced operator is 
written in terms of a distinct set of generalized coordinates corresponding to the local RGBs associated with 
each reduced-order model (ROM). To address this issue, approaches based on congruent transformations 
are proposed in [iHlEniED when the underlying HDMs are defined on a common mesh. These approaches 
cannot, however, be applied when each HDM is defined on a different mesh. The present paper introduces 
a novel approach, also based on congruent transformations, that addresses the challenge associated with 
arbitrary underlying meshes. 

Special care is also given in this paper to the interpolation step of the proposed procedure. The preserva¬ 
tion of properties associated with the linear reduced operators can indeed be enforced by interpolating these 
operators on appropriate matrix manifolds [iniHaHHHUED- In that case, after appropriately mapping the 
reduced operators, interpolation can be carried out in the tangent space to the matrix manifold. As such, as 
long as the interpolation procedure preserves the tangent space, the interpolated quantity will also belong 
to the tangent space and can be mapped back to the manifold, leading to an interpolated reduced quantity 
that preserves its properties. 

This paper is organized as follows. The problem of interest and the ROM database approach are formu¬ 
lated in Section!^ The issue of consistency of reduced-order operators and its enforcement in the case of 
common and arbitrary underlying meshes are then investigated in Section The interpolation of the result¬ 
ing consistent database of ROMs is developed in Section]^ Special care is given to the sampling, storage and 
exploitation of the database. The proposed approach is applied in Section to the model reduction of two 
parameterized systems. The first one is the analysis of a parameterized acoustic scattering system defined on 
arbitrary underlying meshes. In that case, the database approach is applied to the online solution of inverse 
problems. The second problem is the real-time flutter analysis of an aeroelastic system for flight conditions 
ranging from the subsonic to supersonic regimes. It is shown that the proposed approach successfully enables 
real-time predictions on a mobile device. Finally, conclusions are given in Section 


2. Problem formulation and solution approach 

In this paper, linear-time invariant parametric (LTIP) systems of one of the following two forms are 
considered 

1. First-order LTIP systems of the form 

^ A(^)w(i) + B(/i)u(t) (1) 

and their formulation in the frequency domain 

(ju;E(/[x) - A{/j.))w{uj) = (2) 

The high-dimensional state vector is w G = — 1, t > 0 denotes time and cj > 0 frequency. E 

and A are square high-dimensional matrices of dimension N. u G denotes the input variable of 
dimension Ni N and B G ^\\ operators depend on a vector of parameters fi <Z 

For both formulations, an output quantity of interest y G is defined as 

y = G(/[x)w + H(/[x)u, (3) 

with Ao < A and G G and H G 
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2. Second-order LTIP systems of the form 

K(/i)w(i) = B(^)u(t) (4) 

and their equivalent formulation in the frequency domain are considered 

+ jujC{/j,) + K{/j,))w{uj) = (5) 

M, C and K are also parameter-dependent square linear operators of dimension N. For both formu¬ 
lations, an output quantity of interest y can be defined as in (|^. 

The problem of interest is the fast computation of the output y for a given value of the parameters fj, G 
V. More specifically, the computational cost associated with these computations should not scale with N 
anymore. 

To address this problem, an approach based on a database of linear projection-based ROMs is considered. 
This approach proceeds in two steps. 

1. In the first offline step, A^db sample parameter values G P C are selected and ROMs are 

constructed for each parameter value by defining right and left reduced-order bases (ROBs) V(/L6) G 
I^ATx/c W{fjL) G k N and approximating the state w as w « V(/LA)q where q G is 

either solution of a reduced LTIP system in the time domain or frequency domain, as detailed below. 
The output equation in terms of the reduced variable q is then 

Yr = Gr{lji)q + H(/L6)u (6) 


with Grifi) = G{fi)v{^) G 

(a) For first-order LTIP systems, the reduced equations in the time-domain are 

Er(M)^(i) = Ar{n)ci{t) + Br{lJ.)u{t) (7) 

and in the frequency domain 

(juiEriiJ.) - Ar{n))q_{u!) = Br{n)u{u). (8) 

where E^(/Lt) = W{n)'^E{iJ.)V{n) e Ari/u.) = W(M)'^A(/Li)V(/Lt) e and Br{n) = 

(b) For second-order LTIP systems, the time-domain reduced equations are 

W + Kr(M)q(i) = Br{fJ,)u{t) (9) 

where Mr{n) = W(/Lt)^M(/u)V(/x), C^(/Lt) = W(/Lt)^C(/x)V(/u) and K^(/u) = W(/u)'^K(/u)V(/Lt) 
are square reduced operators of dimension k. The reduced equations in the frequency domain are 

+ K^(/L4))q(u;) = (10) 

There are several model reduction techniques that can be applied to construct the ROBs W(/L6) and 

V(/L6) for LTI systems. Among those, the most popular are proper orthogonal decomposition (POD) O 
[22] , balanced truncation [1] and moment matching [23l [24| . 

Once the reduced operators are computed for the sampled values of the parameters, these reduced 
matrices are stored in a database of the form 

VB = {Mi, (E,(Mi), A,(Mi),B,(Mi), G,(Mi), H(Mi))}tT (H) 
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for first-order systems and 

for second-order systems. 

In this work, the left and right RGBs are assumed to have orthonormal columns with respect to a 
common symmetric positive definite matrix Al, that is WAiW{fi) = and V(/l6)^A1V(/l6) = 
I/c. This property can be easily enforced a posteriori by applying a Gram-Schmidt orthogonalization 
procedure to the columns of non-orthonormal RGBs or directly in the RGB construction procedure [25] . 

2. In the online phase, for a given value e V of the parameters, reduced operators are constructed by 
interpolation of the elements of the database VB. 

Two technical issues however arise in the online interpolation step associated with the proposed approach: 

1. The reduced quantities are not defined in the same system of reduced coordinates. Indeed, for each 
parameter /L6^, the system of reduced coordinates is defined by the local RGBs V(/L6j and W(/L6j. 
A naive interpolation of the reduced operators may result in interpolation quantities that are not 
consistent with each other. A comprehensive approach to address this issue is proposed in Section 

2. The linear operators stored in the database VB may have properties that should be preserved by 
interpolation. An approach relying on interpolation on a matrix manifold has been proposed in HZIET] 
to preserve these properties. The approach is briefly recalled and extended in Section]^ 


3. Consistency between reduced-order models 


3.1. Concept 

As underlined in the previous section and in [21], the fact that local reduced operators are defined in 
different sets of generalized coordinates prevents their direct interpolation. In this paper, two approaches 
are presented to address this issue. Both rely on a congruence transformation of the reduced operators to 
enforce consistency. They recognize the fact that the choice of local RGBs is not unique. Indeed, for a 
given right RGB V(/L4), any RGB of the form V(/L4)Q with Q^Q = I/, defines an equally valid coordinate 
representation for the same RGM [21] with Al-orthogonal columns. Similarly, for the left RGB W(/L6), any 
left RGB of the form W(/L6)Z with Z^Z = I/, defines an equally valid basis. 

In turn, for a given first-order LTI RGM IZ = (E^, A^, B^, G^, H), an equivalence class of RGMs under 
left and right multiplications by orthogonal matrices Z and Q can be defined as 

C(7^) = {(Z^E^Q, Z^A^Q, Z^B^, G^Q, H) such that Z^Z = I, Q^Q = l} . (13) 

Similarly, for a second-order LTI RGM IZ = (M^, C^, K^, B^, G^, H), the equivalence class of RGMs under 
left and right multiplications by orthogonal matrices Z and Q is 

C(7^) = {(Z^M^Q, Z^C^Q, Z^K^Q, Z^B^, G^Q, H) such that Z^Z = I, Q^Q = l} . (14) 


Both approaches proposed in this paper will rely on a preprocessing step for which optimal transforma¬ 
tions and Z*(/L4j, i = 1, • • • , Adb are applied to the Adb RGMs stored in the database VB. 

The first approach, originally introduced in m and described in Section |3.2| is applicable whenever the 
underlying HDMs are defined on the same reference mesh. The location of the mesh nodes is potentially 
parameter dependent but the topology of the mesh is common to all parameter values. This requirement is 
relaxed in Section [33] where a novel approach is introduced to enforce consistency in the case of arbitrary 
meshes. In that case, each HDM can have a different number of dofs. 

A word of caution should however be formulated regarding interpolation of reduced linear operators. 
There are cases for which consistency cannot be enforced. Consider for instance the case of two configurations 


4 




fjii and /J .2 associated with a common mesh and for which the right RGBs and ^{/J' 2 ) orthogonal 

to each other. In that case, the subspaces respectively defined by the RGBs are orthogonal as well and no 
transformation of the form {V can define a consistent set of reduced coordinates. The degree 
of consistency between two RGBs will be quantitatively defined in the case of common underlying meshes 
in Section 3.2 and a truncation procedure introduced to further enforce consistency. 


3.2. Enforcement in the case of a common underlying mesh 

Consistency can be enforced in the case of a common underlying mesh by solving a series of Procrustes 
problems [21]. More specifically, given two local right reduced bases and \j = V(/L6^ ), \j can 

be written in terms of as 

V,=VRy+Ti,- (15) 

where Tij is the component of Vj that is Al-orthogonal to i.e. = 0. 

The subspace angles between the RGBs and Vj define a measure of the maximum achievable consis¬ 
tency. The subspace angles can be computed by the following three-step procedure: 

1. Form VfMVj = Kij. 

2. Compute a singular value decomposition where X = [xi,..., x/.], Y = [yi,..., yj^] and 

X = diag(cri, • • • ,cr/e). 

3. Compute the subspace angles as Oi = arccos(cr^), ^ = 1, • • • , /c. The canonical vectors associated with 
each angle Oi are (V^x^, V^y^). Note that the angles are ordered increasingly as 0 < 6>i < • • • < 6^/^ < f. 

A subspace angle Oi that is equal to zero reflects perfect consistency between the associated vectors 
V^x^ and V^y^. In general, angles that are greater than a threshold 6>max = f may define cases for which 
consistency cannot be achieved. Gne option to address this issue that is discussed below is to truncate the 
directions associated with those large angles. An alternate option is to refine the database until smaller 
subspace angles are achieved. 

For a given database PS, optimal transformations can be computed by fixing one of the 

RGBs (say Q(/L6^^) = 1^) and computing the transformations as the minimizers of the following series of 
Procrustes problems: 


Qifii) = argmin||ViS - i = !,■■■ ,Ndb, 

s.t. S^S = Ifc, 

where ||N||jv 4 = The optimal transformation Q(Atj) can be determined analytically from the 

SVD of RiiQ as 

Q(Ai,) = XY^. (17) 

As stated above, truncation of the RGBs can be used to enforce consistency. For each RGB V^, i = 
1, • • • , Ydb, the maximum index £i for which the subspace angles with are smaller than ^max can be 
determined and all RGM truncated to the index 


L = min ii. (18) 

Truncating the RGMs will result in a loss of accuracy of each RGM when compared to the underlying HDM. 
However, this truncation step will improve accuracy after interpolation of the elements of the database as 
these will be more consistent. 
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3.3. Enforcement in the case of arbitrary underlying meshes 

In the case of arbitrary meshes, subspace angles cannot be defined as the underlying HDM spaces may be 
of different, parameter-dependent dimensions N{/j.^). To address this issue, a heuristic procedure is developed 
in this section to enforce consistency for that specific scenario. As in the case of the Procrustes problem, one 
of the ROMs Tli^ is selected to define a reference configuration. Then, for each ROM 7^^, i = 1, • • • , A^db 
in the database, a transformed ROM 7Z^ is determined as the minimizer of a measure of distance of the 
reference ROM IZi^ to the equivalence class C{lZi) as 

TVl = argmin ( 19 ) 

necijZi) 

The minimization problem is schematically depicted in Fig. 



Figure 1 Schematic description of the minimization problem for consistency enforcement in the case of arbitrary 
meshes 

The measure of distance is defined as follows: 

1. For a first-order system, the distance between IZ = (E^, A^, B^, G^, H) and VJ = (E^, A^, B^, G^, H') 
is defined as the normalized expression 

Di,{n,n') = e\\Er - KWl + «I|A. - KWl + /3||B, - B'Jl + 7||g, - g;ii|. + 7/||H - H'lll, (20) 


where 


^ 1 1 

" IIE0|iy “ i|A0|iy ^ 


2’^ "B0|||’ ^ l|G0||y ||H0||2, 


( 21 ) 


are normalization constants based on the reduced operators in IZi^ = (E^, A^, Bj:, G^, H^). 


2. For a second-order system, the distance between IZ = (M^, C^, K^, B^, G^., H) and IZ' = (M^, C^, K^, B^, G^, H') 
is defined as 


r' l|2 


Ao(^,7^') = Ml|M,-M;|||.+C||C,-C;|||. + «||K,-K;|||+/?||B,-B;|||.+7||G,-G;|||.+^?||H-H'||^ 

( 22 ) 

where 

1.1 1 


h = 


■r IlF 






K, = 


r IlF 


(23) 


are 


normalization constants based on the reduced operators in TZi^ = (M^, C^, K^, B^, G^, H^). 
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In practice, since the class C(lZi) is parameterized by two transformation matrices Q and Z, the optimiza¬ 
tion problem ( p^ can be explicitly written for a first-order system in terms of IZi = (E^^, H^), 

IZiQ and (Q, Z) as 


mm e 

Q,Z 


|Z^EhQ-E°| 


%+a\\ 


Z^ArQ - A°|||. + PWZ^Br - B« 


II 


G,Q-G0||i. + 7?||H-H0| 


s.t. Q Q = I/c, Z^Z = I/c. 


(24) 


A similar expression holds for second-order systems. The rest of this section will focus on first-order systems 
but the analysis directly carries over to second-order systems as well. 

Remark. In the case of Galerkin projection, W(/L4j = V(/L6j and Q = Z. Then (19) simplifies to 


min ellQ^EnQ - E0|||. + a||Q^A,Q - A°||| + ^||Q^B, - B°||| + 7||G,Q - G0||1. + 7?||H - H0|||. 

v4 

S.t. Q^Q = I/c. 

(25) 


Problem (24) is equivalent to the maximization problem 


max (eZ^ErtQ, E°) + (aZ^A^Q, A°) + {/3Br (B°)^ , Z) + ( 7 (G°)^G^, Q) 

Q ,Z 

s.t. Q^Q = Ife, Z^Z = Ik, 


(26) 


where 

(M,N) =tr(M^N), (27) 

The problem of maximizing the first term in Eq. ( |26| ) has been studied in the literature [261127] in the 
case Q = Z (Galerkin projection). This first term defines a correlation criterion between the matrices 
F/ri and E^. A solution to that problem developed in Ref. [26], consists of defining an iterative algorithm 
whose fixed points are the critical points of the maximization problem. This approach is here extended to 
the optimization problem of interest for enforcing consistency between ROM operators. Both Galerkin and 
Petrov-Galerkin projections are considered as follows. 


Galerkin projection. 


In this case, the functional in (26) is 


Jg{Q) = (eQ^EnQ, E°) + (aQ^AnQ, A°) + (/?B, (B°)^ , Q) + (7(G°)^G,, Q). 


(28) 


Adapting the algorithm developed in [26j to the present case, the iterative algorithm is based on an 
affine map defined as 


M,,g(Q) = e (EhQ (EO)^ + E^,QEo) + a (a^Q (A^)^ + A^,QA°) + sQ + pB^i (B^)^ + 

(29) 

where s is a fixed real parameter chosen such that s > 5min,G with 


s„,i„,G = 2e ||Eh ||2 ||E °||2 + 2a ||Ah ||2 HA^H^ + /3Bh (B°)^ + 7 G^,G° 


(30) 


Defining the parameter s is necessary to ensure that the fixed points of the proposed iterative algorithm 
are exactly the critical points of the maximization problem (see Theorem 1). 

The proposed procedure then proceeds by iteratively solving the maximization problem 


Qi+i=arg max (S, Ms,g(S 2 )), j = 0, ■ • • • 

S^S=Ifc 


(31) 
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The solution to this problem is established in the following lemma, leading to the proposed iterative 
procedure presented in Algorithmic 

Lemma. Let the singular value decomposition of 'WLs,G{Qj) be defined as 

M,,G(Qi) = USV^. (32) 

Then 

k 

^ (Q,M^(Qj)) = ^cr<!(M^,G(Qi)), (33) 

Q Q=Ifc 

where {cr^(Ms^G(Qj))}t=i fbe set of singular values of Ms^G(Qj)- The solution to the maximization 
problem is unique and equal to UV^. 

A proof is offered in [26] in a more general setting. 


Algorithm 1 Fixed-point procedure in the case of Galerkin projection 

1: Compute s>2e ||E^i ||2 + 2a \\Ari \\2 ||A ^||2 + 

2 : Compute F = (B^)'^ + tG^^GO. 

3: Choose an orthogonal initial matrix Qo G 
4: for jf = 0, • • • do 
5: Compute the map 

M,,G(Qi) = e (EhQ, (E°)^ + E^,Q,E°) + a (AhQ,- (A°)^ + A^,Q,.A°) + sQ,- + F 
6: Compute its SVD 

Uj+iSj+iVj_,_i = ]V[5^G(Qi) 


pBri (B0)^ + 7G^,G0 


7: Qi+1 — 

8 : end for 


In order to show that the fixed points of the recursive algorithm are the critical points of J'g^ one needs 
to characterize these critical points. This is done in the following theorem. 

Theorem 1. The critical points Q"* *" of Jg are orthogonal matrices satisfying the identity 

Q*S = e (EhQ* (E°)^ + E^iQ*E°) +a (a^Q* (A°)^ + A^.Q^A^) +/3Bh + (34) 

where S is a symmetric matrix. 

A proof is offered in Appendix 1. 

Theorem 2. The set of the fixed points of the recursive Algorithm [C is exactly the set of the critical 
points of J'g- 

A proof of the theorem is presented in Appendix 2. 

• Petrov-Calerkin projection. 

Defining the functional 

Jpg(Q,Z) = (eZ^ErtQ,E°) + {aZ"^AriQ,A°) + (/3B, (B°)^,Z) + (7(G°)^G,,Q), (35) 
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the iterative algorithm is now based on the block-affine map defined as 
Ms^pg(Q, Z) = 


M^pg(Q,Z) 

0 



0 

M,^pg(Q,Z) J 


■ 0 Eri ' 

■ Q o' 

1- 

Sh 

O 

& 

O 

1 _ 

.Ki 0 

0 Z 

. (E“ 

f 0 J 


a 


0 A.j'j 

■ Q o' 

0 

1 

O 

_Al 0 _ 

0 Z 

O 

_1 

0 


' Q o' 

_ 1 _ 

■ (BO)^ 0 

0 Z 

1 

0 7G^^G° _ 


where s is chosen such that s > Smin pg with 


«Smin,PG = e ||E, 


e; 


01 


^*112 ||■^r||2 


+ G II A, 


■rz\\2 


r\\2 


+ max 


(^||Brt (BO)' 


( 36 ) 




.) ■ 


Similarly as in the case of Galerkin projection, a fixed point procedure is defined in Algorithmic The 


Algorithm 2 Fixed-point procedure in the case of Petrov-Galerkin projection 


1 : 

2 : 

3: 

4: 

5: 


Compute s > e ||E^i ||2 IjE^H^ + o ||A^i ||2 \\^^2 


max 





Compute Fb = pBrt (B^)^ and Fg = tG^^GO. 

Choose orthogonal initial matrices Qo G and Zq G 

for J = 0, • • • do 

Compute the maps 


(b;)’'i|j,7||gsg; 


,)■ 


6 : 


^^PG(Qj’ ~ eE^^Zj (E^) -f- oA^^Zj (A^) + sQj + 

and 

= eE^,Q,- (eO)^ + aAj,Qj (A^)^ + sZj + Fb 


Compute their SVDs 


and 




7: 

8: Z,+i=Uf+iVf7i 

9: end for 


following theorems, whose proofs follow closely the ones of Theorems 1 and 2, establish the fact that 
the fixed point procedure in Algorithm [C can be used to find critical points of Jpg- 

Theorem 3 The critical points (Q"*", Z*) of Jpg are orthogonal matrices satisfying the identities 

Q*Sq = cEhZ* (eO)^ + aAriZ* (A^)^ + 

Z*Sz = eE^,Q* (E°)^ + aAj,Q* (A°)^ + / 3 Bh (B°)^ , 
where Sg and are symmetric matrices. 

Theorem 4. The set of the fixed points of the recursive Algorithm |C is exactly the set of the critical 
points of Jpg- 
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Remark. For the case of arbitrary meshes, assessing consistency is a more difficult task as two conflicting 
factors intervene in the distance measure Di^: (1) the inconsistency arising from a choice of two distinct sets 
of coordinates and (2) the inherent variation of the ROM operators due to parameter changes. 

3.4- Consistent set of redueed-order models 

After the computation of the optimal transformation operators the reduced op¬ 

erators in the database VB are transformed accordingly as 

VB = (E*(Mi), A*M, B*M, G*M, 

(39) 

Similar expressions hold for the case of second-order systems. 

4. Interpolation in a database of ROMs on matrix manifolds 

4-1. Interpolation 

As indicated in Section the interpolation of the linear operators stored in the database VB should 
often preserve properties of the operators such as symmetry, positivity, orthogonality or non-singularity. An 
approach to preserve these properties was first presented in m- It is based on the interpolation on the 
tangent space of the appropriate manifold and was applied to the case of interpolation of reduced operators 

in [IIlIIiET]. 

The algorithm proceeds for each of the elements of the database in four steps as follows: 

1. An identification of the manifold the reduced matrices belong to 

2. A mapping (logarithmic map) of all the database reduced matrices to the tangent space of the manifold 
at one of the database points 

3. An interpolation of the mapped quantities in the tangent space at the target parameter fj, 

4. A mapping (exponential map) of the interpolated quantity back to the manifold leading to a reduced 
operator at the target parameter fj, 

More details on the interpolation algorithm as well as the formulas for computing the mapping are provided 

in EH HI]. 

In practice there may be several choices for an interpolation procedure on matrix manifolds as underlined 
by the following two cases. 

• In [19] , the authors develop a heuristic technique for interpolating non-singular matrices either on the 
manifold on non-singular matrices or square matrices. The heuristic proceeds by selecting the manifold 
for which a nonlinearity criterion is the smallest. This heuristic is applied in Section [53] 

• An alternative to interpolating symmetric positive definite matrices on the tangent space to that 
manifold is to use the Choleski factorization. This novel approach is described in Algorithm and 
applied in Section [5Tj This approach avoids selecting one of the database points and interpolating in 
its associated tangent space. It preserves the SPD properties of the matrices as long as the interpolated 
quantity on the diagonal of the Choleski factor are all strictly non zero. 

There is no restriction on the interpolation technique in the tangent space to the matrix manifold of 
interest as long as the it leads to an interpolated quantity that preserves the tangent space EH- [28] . 
the author identifies an interpolation technique that does not preserve that property. When the database 
parameters belong to a lattice of points, spline or polynomial interpolation can be used in the tangent 
space. When the dimension of the parameter domain is large, however, interpolating from a lattice of 
points is not an option and instead, interpolation based on radial basis functions or Kriging can be used 

instead uni nails]- 
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Algorithm 3 Interpolation of SPD matrices by Choleski factorization 

1: for i = 1, • • • , Adb do 
2: Compute the Choleski factorization 

3: end for 

4: Interpolate the Choleski factors leading to an interpolated factor S* 

5: Compute the interpolated matrix as 


4 . 2 . Sampling 

The selection of sample points is an important step that influences the accuracy of the resulting 

interpolation approach. A poor choice of sample points will result in large errors of the proposed procedure in 
some regions of the parameter domain V. A priori sampling techniques such as factorial and latin hypercube 
sampling can be used to provide a uniform coverage of V. Alternatively, greedy techniques that iteratively 
sample the regions of the parameter space associated with the largest ROM error can provide a selection of 
the samples that is more suited for the problem of interest. Such greedy techniques have been introduced 
in the context of model reduction in general in [9l [30l |3T1 [131 HH El] and for interpolation of LTIP ROM 
systems in particular in [29]. A priori sampling will be used in the application of Section 5.2 and a greedy 
sampling approach developed for the inverse problem application of Section [5Tj 


4 . 3 . Storage and exploitation 

In practice the reduced operators are stored after their congruence transformation in one database VB 
or several sub-databases of the form 

Ns Ns 

VB=[jVBs=[j (40) 

S = 1 S = 1 

for first-order systems and 

Ns Ns 

vB=[jvB,= \J (4i) 

S = 1 S = 1 

for second-order systems. 

Storing the database is inexpensive as it only involves reduced operators. In practice, a database with 
^DB = ^DB,s contains + k{Ni + No) + NiNo) matrix entries for first-order systems of 

the forniQ and Ndb{N^ + 3k‘^ + k{Ni + No) + N^No) entries for second-order systems of the form Q. 
The parameter domain V is in practice subdivided in Ns non-overlaping subdomains 

Ns 

V=[jVs (42) 

S = 1 

and each subdomain Vg is associated with the sub-database VBs- Then, in the online phase, for a new value 
/i G P of the parameters, the sub-database VBsq it belongs to is readily identified and reduced-operators 
computed for using the ROMs stored in VBsq • 


5. Applications and performance assessment 

The proposed approaches are here applied to two challenging physical applications: the acoustic inverse 
obstacle problem and the aeroelastic flutter problem. 
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5.1. Acoustic scattering analysis 

The acoustic inverse obstacle problem considered here consists in determining the shape of an obstacle or 
a part of this shape from the knowledge of some scattered far-field patterns, assuming certain characteristics 
of the surface of the obstacle. It is well-known [33] that such an inverse problem is non-linear and often quite 
ill-posed, making its numerical solution challenging. 

To illustrate the ROM database framework proposed above, a parameter identification problem is consid¬ 
ered, where the shape of the obstacle is assumed to be known a-priori, but the vector of parameters ijl G 
of the shape needs to be identified from the measured far-held pattern. This class of problems is a subset of 
a general acoustic inverse obstacle problem. 

In order to describe the considered problem more accurately, the corresponding direct acoustic scattering 
problem is recalled hrst. The scattering of time-harmonic acoustic waves by an impenetrable obstacle with 
the boundary S embedded in an inhnite homogeneous huid medium C can be formulated as the 
following exterior boundary value problem for the unknown acoustic pressure held w in the huid 

Are + = 0 in 

(a + b-^ (w + w*”") = 0 on 

ctzi (dw . \ 

hm r 2 —- jnw ] = 0 , 

r^oo y OT J 




(43) 


where the incident wave is given by 


W^nc = gi/^d-x 


the unit vector d indicates the direction of the incident plane wave, and either a 7 ^ 0 or 6 7 ^ 0 . Sound-hard, 


sound-soft, or impedance boundary conditions can all be represented by the second equation of (43). In the 


example below, the sound-hard scattering problem is used, leading to the choice of the Neumann boundary 


condition (a = 0 and 6 = 1). The third equation in (43) is the Sommerfeld radiation condition. It ensures, 
in the physical sense, that all waves are outgoing and, mathematically, that the direct scattering problem is 
well-posed for any wavenumber n = cj/c, where uj is the angular frequency of the harmonic oscillations and 
c is the speed of sound in the fluid. 

In order to discretize the direct scattering problem (43), the finite element method is considered. The 
infinite domain is first truncated, and a perfectly matched layer [34] near the exterior boundary is used to 
simulate the effect of the Sommerfeld condition. This converts the boundary value problem (43) into that 
of solving the algebraic system of linear equations 


(K(|i) - K2]V[(|i)) w(K,|i) = f(K,/Li) (44) 

for the unknown degrees of freedom w e . This system is a second-order LTI system of the form (|^ . Here, 
K corresponds to the finite element discretization of the Laplace operator, and M is a mass-type matrix. For 
an interior problem associated with the Helmholtz equation, M is real, symmetric positive-definite, and K 
is real, symmetric non-negative. When the perfectly matched layer is used, the matrices K and M become 
complex and non-Hermitian. The source vector f arises from the discretization of the sound-hard boundary 
condition. 

The far field pattern characterizes the asymptotic behavior of the acoustic scattered field far away from 
the obstacle. In two dimensions, it admits the following integral representation [33] 


Wr 


>(x) = 


p3 4 


r J f ^(y) + ^(y) 1 ® 


{Sirnr Jr Kdi' 


&S\ 


(45) 


where is the unit circle. After computing the finite-element solution, the integral (45) can be evaluated 


by integrating over a suitable curve T (often the boundary S) in the computational domain. The integral in 
(45) evaluated at Nq locations xi, • • • , of the circle can then be in practice represented by the action 
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of a matrix on the solution vector y(/^,/L4) = fi) with G G Practically, the following 

logarithmic scale quantity is usually plotted 

5(x) = 101ogio(27r|«;oo(x)p). 

This quantity can also be computed for each entry of the output vector y(/^, /l^) as 

s(k,m) = s(y(K,Ai)) . (46) 

The inverse problem considered here consists of identifying parameters fi = G C of a 

two-dimensional mockup submarine, characterized by its length and the position of its tower from 
given far-field data for several frequencies (wavenumbers) The inverse problem can be written 

as 


min - Sm{i^i 




S.t. (K(^) - KfM(/i)) w(Ki, ^l) = f (Ki, ^l), 

y(Ki, ^l) = G(Ki, ^l)w{Ki, ^l), i = !,■■■ , N^, 


(47) 


where a Tikhonov regularization term has been added and l3 are appropriate positive weights. 

Since the solution of the discrete direct problem (44) for each different wavenumber requires a costly re¬ 


factorization of the matrix on the left-hand side, to enable the efficient computations for many wavenumbers, 
a reduced-order model is built using a derivative-based Galerkin projection (DGP) framework [24] for a given 
value fi of the shape parameters. In this method, based on moment-matching, derivatives of the 

solution with respect to the wavenumber k, are first computed by solving with recursively constructed 
right-hand sides at interpolating wavenumbers. Then, these derivatives are orthogonalized to achieve 

numerical robustness, and used to form a subspace of dimension k = for Galerkin projection, 

leading to reduced matrices K^(/L6), M^(/L6), and 

Figure [^depicts the triangulated computational domain (left) with the elements in the PML layer shown 
in cyan; the real part of the solution for /i: = 20, = 1, and = 0.2 is shown on the right. For different 

values of the shape parameters, the computational domain is remeshed. Isoparametric cubic Finite Elements 
are used. All computations are done using Mat lab. 

The solution procedure described in detail below builds a database of frequency-sweep ROMs offline by 
sampling the shape parameter space adaptively to ensure accuracy. The database of the ROMs is then used 
to efficiently solve the reduced inverse problem. For a given value of the parameters, = 8 (including 

the 0-th derivative) are computed for the = 2 frequencies /i: G {10, 20}, leading to a ROM of dimension 

k = 16. In the present case, there are Nq = 360 outputs equidistributed on the sphere S^. 

In a first set of numerical experiments, the effect of consistency on predictions based on a database of 
ROMs is illustrated. For that purpose, a small database of A^db = 4 ROMs is first created for the values 
of the parameters indicated in Table where the number of dofs for each underlying HDM is also reported. 
One can observe that each HDM has a different number of dofs. 


The consistency enforcement procedure for arbitrary meshes developed in Section 3.3 is first applied 
to transform the reduced operators prior to their interpolation. In a second case, the operators are not 
transformed. In all cases, the reduced operators are interpolated on their appropriate manifold at = 
[0.9625,0.1125] by bilinear interpolation as follows: 


• The operators Re(M^(')), Im(M^(')) and Re(K^(')) are SPD matrices and as such as interpolated on 
their appropriate manifold using the Gholeski decomposition-based approach proposed in Section |4.1[ 

• The operator Im(K^(-)) is interpolated on the manifold of symmetric matrices. 

• The operators {fr(/^ 2 , •), Gr(/^i, are interpolated on the manifold of rectangular complex matrices. 
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Figure 2 Computational mesh and the real part of the solution for = 20 


VB point 

/^i 

^2 

7V(m) 

1 

0.95 

0.1 

41,235 

2 

0.975 

0.1 

40,965 

3 

0.95 

0.125 

41,424 

4 

0.975 

0.125 

40,929 


Table 1 Database of ROMs and associated underlying HDM number of dofs 


Figure [^reports the far-held pattern for n = lA dX = [0.9625,0.1125] computed by the HDM and the 
two ROM interpolation approaches. One can observe the importance of consistency as the inconsistent ROM 
database leads to poor predictions whereas the consistent ROM database predictions very closely follow the 
HDM results. 

Next an adaptive approach for constructing the ROM database is developed. The approach proceeds by 
comparing the predictions associated with the ROM database and HDM at the center of each hypercube of 
the database and rehning that hypercube if the error is above a given threshold. The error can be associated 
with an output of interest, such as the far-held pattern, leading to a goal-oriented approach. In the following, 
the error will be based on the accuracy of the shape returned by solving the reduced inverse problem by 
simulated annealing. 

i=i (48) 

S.t. (Kr(/i) - K^Mrifl)) q(Kj, fl) = triKi, n), 

Yr 5 M) ~ : M) : i = 1, * * * , N ^. 

The error measure is then dehned as 


Error = 


fl- 


Mmax Mmin 


CX3 


(49) 


where denotes the target shapes, A^rnin and A^rnax ^^le lower and upper bounds for each shape parameters, 
respectively and the ratio in Error is computed entry-by-entry. 


14 




































Figure 3 Comparison of the far-field pattern predictions obtained by interpolation of inconsistent and consistent 
databases. 


Figure reports the refined database for an error tolerance of 5% and the parameter domain V = 
[0.9,1] X [0.1, 0.2]. A^db = 21 points are sampled in the domain. The training errors obtained at each 
iteration refinement of the procedure are reported in Figure One can observe that after the second 
refinement, all errors are below the error threshold of 5%. 

To validate the accuracy of the ROM database, 289 target shape parameters are selected in V and the 
reduced inverse problems solved for each of them. The distribution of corresponding errors are reported in 
Figures!^ andOne can observe that all errors are below the error threshold, confirming the validity of the 
training procedure. 

Finally, the CPU timings associated with the solution of a given inverse problem are compared for the 
HDM and ROM database strategy. The CPU timings are reported in Table One can observe that an 
impressive speedup of 270 is obtained thanks to the database strategy. For a given function call in the 
optimization problem, the speedup is equal to 207. 


Approach 

Error 

Number of 
function calls 

Online optimization 
Wall time 

Speedup 

HDM 

8 X 10“^ 

1530 

1 h 30 min 

1 

ROM database 

0.02 

1176 

20.1 s 

270 


Table 2 Wall times associated with the solution of the inverse problem with the HDM and database of ROMs 


5.2. Flutter analysis 

The aeroelastic analysis of a wing-store configuration flying in the subsonic, transonic and supersonic 
regimes is considered. Some properties of that system were originally studied in [35] [36]. Among those, 
the hydroelastic effects inside the fuel tank modify the structural properties of the wing-store configuration, 
thereby affecting the flutter characteristic of the system. This justifies parameterizing the aeroelastic system 
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Pre-computed database VB of reduced-order models selected by the adaptive sampling procedure. 
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Error 



Error 

Figure 6 Distribution of prediction errors for the solution of 289 inverse problems. 



0.1 0.9 /ii 

Figure 7 Prediction errors for the solution of 289 inverse problems. 
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(a) CFD surface grid (b) FE structural model 

Figure 8 High-dimensional aeroelastic model of a wing-store configuration. 


by the fuel fill level / inside the tank. Furthermore, the aeroelastic properties of the wing-store system 
depend heavily on the aerodynamic properties of the configuration. As a result, the system will also be 
parameterized by the free-stream Mach number M^. 

The structural and fluid surface models of the wing-store system are graphically depicted in Figure 
The structural subsystem is a second-order LTIP of the form and is modeled by the Finite Element 
method. For each fill level, a new structural mesh is generated inside the tank for the full physical domain. 
The hydroelastic effects are modeled by an added mass effect [35l |36], resulting in a linear HDM with 
A^(^) = 6, 834 dofs for all values of /. The proposed ROM database approach will enable by-passing the 
re-meshing of the fuel domain every time the fill level is varied. 

The fluid subsystem is modeled by the linearized Euler equations and discretized by the Finite Volume 
method using a second-order accurate linear flux reconstruction and a second-order accurate implicit back¬ 
ward difference time integration scheme. For each operating point (jl = (Moo, /), the nonlinear fluid HDM is 
linearized around a steady-state, resulting in a first-order LTIP system of the form Q with ~ qpQ, 000 
dofs. 

In this work, the operation domain of interest is (Moo, /) G P = [0.6,1.1] x [0,100]. For each operating 
point, the critical values of pressure and velocity at the onset of flutter are sought. Once, these 
quantities are determined, the flutter speed index (FSI) can be computed as 


ycr 

FSI = 




(50) 


where hg is the semi-chord of the wing at its root, cJq, is the first dry torsional mode of the wing-store 
structural system and p is the mass ratio as defined in [3711351 [361. 

The flutter speed indices of the system of interest are computed using the HDM for 26 different free- 
stream Mach number and 5 different fill levels in the domain (Moo, /) ^ 7), resulting in 130 operating points, 
and reported in Figure One can observe the characteristic flutter dip for Moo ~ 0.96. 

The framework developed in this paper is then applied to the problem of interest to interpolate reduced 
aeroelastic operators. In this example, all structural and fluid HDMs are defined on the same mesh and the 


approach of enforcing ROM consistency developed in Section is followed. For that purpose, Vdb =21 
operating points are sampled and their corresponding aeroelastic ROMs constructed and stored in the offline 
phase in a database VB. These points correspond to a lattice (Moo, /) ^ {0-b, 0.75,0.9,0.95,1.0,1.05,1.1} x 
{0,50,100}. For each operating point, a structural ROM of dimension = 4 is constructed by projecting 
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the linear structural HDM onto its first four natural modes. Then a fluid ROM of dimension = 15 is 
constructed by POD using the approach described in [inmiiiiH] and a reduction of the system in descriptor 
form [25] . The FSIs predicted by those 21 aeroelastic ROMs are reported in the right portion of Figurej^ Very 
good agreements can be observed at the database points when compared to their counterparts determined 
by the HDM that are depicted in the left portion of that same figure. 

HDM Database 




Figure 9 Comparison of the high-dimensional model and database reduced-order models flutter speed indices. 

The 21 pre-computed aeroelastic ROMs are then distributed in = 3 sub-databases: the first one 
covers the subsonic and lower transonic flow regime M^o C [0.6, 0.9], the second one the upper transonic 
regime M^o C [0.9,1.0] and the third one the supersonic regime Mo© C [1.0,1.1]. These three databases are 
graphically depicted in Figure [Tol In the online interpolation procedure, in each database, piecewise-linear 
interpolation will be used in the Moo direction and cubic spline interpolation in the / direction. 



# Database 1 
0 Database 2 
V Database 3 


Figure 10 Pre-computed sub-databases {T>13s}l=i of aeroelastic reduced-order models. 
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Before the online phase, the heuristic developed in m and mentioned in Section [41] is applied for choosing 
the manifold onto which to interpolate the reduced fluid operators. For the fluid operators, interpolation 


can indeed be done on the manifold GL of non-singular matrices of size or on Since 

the interpolation procedure of choice involves two points in the direction and three in the / direction at 
a time, the heuristic is applied for six different regions of the parametric space, as indicated in Table As 
reported in Table the manifold GL is chosen in two regions while the manifold chosen 

in four regions. 



Database 1 

Database 2 

Database 3 

Moo e 

[0.6,0.75] 

[0.75,0.9] 

[0.9,0.95] 

[0.95,1] 

[1,1.05] 

[1.05,1.1] 

GL 

0.89 

0.87 

0.93 

0.93 

0.76 

0.86 


0.84 

0.86 

0.98 

0.87 

0.78 

0.82 

Manifold choice 



GL 


GL 



Table 3 Non-linearity indicator in each database for the manifold choice heuristic 



Figure 11 Comparison of eigenvalues of the structural subsystem: HDM □, interpolated ROM o. 


Next, the proposed methodology is applied to interpolate the aeroelastic ROM operators. The properties 
of the structural operators resulting from that interpolation are first analyzed by comparing their respective 
eigenfrequencies to their HDM counterparts. The corresponding results are reported in Figure EH Good 
agreements can be observed, even for fill levels that are not present in the database. Next, the interpolated 
aeroelastic ROM operators are used to predict the onset of flutter in the entire parametric domain (Moo, /) C 
[0.6,1.1] X [0,100]. The predicted FSIs are reported in Figure [l^ Very good qualitative and quantitative 


agreement can be observed for all flight conditions considered. For comparison, response surface estimation 
(RSE) is also applied to predict flutter using the database FSI data reported in the right portion of Figure]^ 
In this case, bicubic spline interpolation is used. When compared with the predictions arising from ROM 
interpolation, the results from RSE are found to be much less accurate, especially near the transonic dip 
and in the supersonic regime. RSE cannot, in particular, predict the ESI behavior for low fill levels at 
supersonic speeds. It is quite remarkable that the method proposed in this paper can capture this complex 
phenomenon with only the ROM database associated with the results shown in the right portion of Eigurej^ 
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Response Surface 


ROM Interpolation 




Figure 12 Comparison of flutter speed indices predicted using (1) response surface estimation, (2) ROM 
interpolation. 


This example underlines the potential for accurate predictions of the proposed method which operates by 
interpolating models and not outputs, as in RSE. The offline and online CPU times associated with the 
prediction of the entire parametric FSI surface each of the four techniques are reported in Table These 
results clearly demonstrate the real-time capability of the proposed approach, as it can accurately predict 
the FSI for 130 configurations in only 31 seconds. 


Approach 

Offline phase 

Online phase 


Wall time (s) 

Number of processors 

Wall time (s) 

Number of processors 

HDM 

- 

- 

9,152,000 

32 

RSE 

28,000 

32 

2 

1 

ROM database 

28,000 

32 

31 

1 


Table 4 


Wall times associated with computing the complete predicted FSI surface reported in Figure and 
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The flutter behavior at supersonic speed for low fuel levels is then studied more in detail by predicting 
the FSI at an additional fuel level for 15% tank fill. The corresponding FSI are reported in Figure 13 In 


order to understand the physical phenomenon at play, the eigenvalues of the interpolated aeroelastic ROM 
matrices are computed for increasing values of free-stream pressure until flutter is reached. By following the 
wet structural modes in the complex plane, one can determine which structural mode flutters, that is which 
one is the first to cross the imaginary axis. Results are reported in Figuresandfor four different values 
of the free-stream mach number and an empty tank. One can observe that for Mo© = 1.075 and = 1.091, 
the first mode is the first to cross the imaginary axis while for = 1.092 and Mo© = 1-1, the third mode 
is the first one to flutter. These results clearly show that a bifurcation phenomenon is at play. Being able 
to perform such analyses demonstrates another clear advantage of the proposed method over RSE. For that 
same fill level, the HDM predicts the same phenomenon, that is a peak of FSI in function of the free-stream 
Mach number M^o between M^o = 1-09 and M^o = 1.092, which is in perfect agreement with the results 
arising from the interpolated ROM. 


Next, the effect of the manifold choice heuristic is studied by considering interpolation on 




and on GL{k^^^), respectively, for the reduced fluid operator at every operating point. The corresponding 
results are reported in Figure for interpolation using Choice 1 and 3, respectively. One can observe that 
the predicted FSI are more accurate when the heuristic is used, especially for subsonic and transonic flight 
conditions. 

Finally, the effect of consistency on the results accuracy is investigated. Inconsistent interpolation for 
both the structural and fluid subsystems and for only one of those two subsystems is performed and the 
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FSI 


HDM 


Response Surface 



0.4 

0.35 

0.3 

0.25 

0.2 

50 


1.1 



ROM Interpolation 



1.06 

Fill Level % 0 1-05 Mach Number 


Figure 13 Comparison of high-dimensional model and predicted flutter speed indices between Moo = 1-05 and 
Moo = 1.1 for low fuel fill levels using (1) response surface estimation, (2) ROM interpolation 
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Figure 14 


Aeroelastic matrix eigenvalues loci at 0% fill level for various free-stream Mach numbers. 


Moo = 1-075 






Moo = 1-091 




-1 -0.5 0 0.5 1 

Real Part 


Figure 15 Aeroelastic matrix eigenvalues loci at 0% fill level for various free-stream Mach numbers (zoom 
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ROM Interpolation 


ROM Interpolation 




Figure 16 Predicted flutter speed indices using ROM interpolation without the manifold choice heuristic: 
interpolation of the reduced fluid operator on (left) and on GL{k^^^) (right). 


ROM Interpolation 


ROM Interpolation 


CO 



0 . 4 ^ 

0.35 

0.3 

0.25 

0.2 


Fill Level % 


Mach Number 



ROM Interpolation 



Fill Level % 0 0-6 Mach Number 


Figure 17 Predicted flutter speed indices using ROM interpolation without consistency enforcement: for the 
structural operators only (top left), for the fluid operators only (top right), for all operators (bottom). 
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corresponding FSI results reported in Figure The reader can observe the crucial effect of consistency as 
none of the interpolation of inconsistent ROMs leads to accurate predicted FSL 

The interpolated ROM can also be used to predict the displacement at a given location of the wing- 
tank system. Here wing tip displacements time histories are predicted using the interpolated ROM at the 
transonic tip, that is M^o = 0.97 for an empty tank. The corresponding results are reported in Figure 
for three different cruise altitudes, and compared to predictions using the HDM. Good agreements can be 
observed. 


M = 970, fill level = 99 %, sea level 



M = 970, fill 

I 


g 0.02 
riS 

0 

Q 

a 

^ - 0.02 

) 

-0.04 


S ■ 


-0.06 


level = 99 %, altitude: 10,000 ft 

" aMi a a 


- HDM 

-Interpolated ROM 


0.2 0.3 

Time (s) 


M = 970, fill level = 99 %, altitude: 20,000 ft 



Figure 18 Wing tip displacement at Moo = 0.97 with full tank predicted using the high-dimensional model 
and an interpolated ROM at various altitudes. 


Finally, in order to demonstrate the capability of the proposed method to operate on mobile devices, an 
iPhone application is implemented for the aeroelastic system of interest. A screenshot of the application 
is displayed in Figure The application can operate in the following two modes, based on the database 
of A^db = 21 points considered in this section: (1) In the first mode, for a given value of the fill level /, 
the FSI is compiled for G [0.6,1.1]. This is the mode depicted in Figure 19 (2) In the second mode. 


for a given combination (Moo,/), the smallest aeroelastic damping ratio is computed for the altitude range 
h G [0,40000] ft. 


6. Conclusions 

This work presents a framework for real-time predictions based on a database of linear reduced-order 
models. It is based on the offline pre-computation of reduced-order models and their online interpolation 
at unsampled values of the parameters. A pre-processing step is first established to enforce the consistency 
of the set of generalized coordinates each reduced operator is defined by. The present paper presents such 
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Figure 19 Screenshot of the iPhone application depicting the FSI for fill level / = 43.77%. 


a step for both systems defined on common and arbitrary underlying meshes. The operators are then 
interpolated on the tangent space to a matrix manifold to enforce properties associated with each operator. 
The framework is then applied to two challenging multi-physics applications, demonstrated its capability to 
lead to real-time and accurate predictions. 
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8. Appendix 1: proof of Theorem 1 

The objective function can be written as 

Jg(Q) = e(Q^EHQ, E°) + ^(Q^AhQ, A°) + (F, Q), 

where F = (BO)^ + 

The Lagrangian of the optimization problem is then 
£(Q, S) = Jg(Q) + Ik - 

= eiQ^EriQ, E°) + aiQ^AriQ, A°) + (F, Q) + ^ is, Ifc - , 


(51) 


(52) 


where G ^ 


cx/c 


is a symmetric matrix of Lagrangian multipliers. Using the following identities 


Vq(M,Q) 

Vq(M,Q^Q) 

Vq(Q^MQ,N) 


M 

Q(]V[ + M'^) 
MQN"^ + M^QN, 
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(53) 

(54) 

(55) 











the gradient of the Lagrangian with respect to Q is obtained as 

Vq£(Q, S) = e (EhQ (eO)^ + E^^QE^) + a (AhQ (A^)^ + A^^QA^) + F - QS (56) 

which leads to the first-order optimality condition 

QS = e (EhQ (E°)^ + E^iQE°) + a (a^Q (A°)^ + A^,QA°) + F, (57) 

together with the constraint Q^Q = Ik and the property that S is symmetric. 


9. Appendix 2: proof of Theorem 2 


The goal of this section is to prove that the set of fixed points of the proposed recursive algorithm is 
equal tojbhe set of critical points of the objective function Jq- ^ ^ ^ ^ 

Let Q denote a fixed point of the recursive method defined in Algorithm 2. Then Q satisfies Q = UV^ 
where 


is a singular value decomposition. Since V is an orthogonal matrix, 


that is 


QS = e (EhQ (E°)^ + E^,QEo) + a (AhQ (A^)^ + A^^QA^) + F, 



J + sQF 

(58) 

y + A^,QA) 

+ sQ + F, 

(59) 

" + A^,QAo) 

+ F, 

(60) 


where S = V5]V^ — is a symmetric matrix. U and V being orthogonal, Q is orthogonal as well and 
therefore meets the requirements of Theorem 2. The set of fixed points of Algorithm 2 is included in the set 
of critical point of Jg> 

Conversely, let Q* be a critical point of Jg> Then, there exists a symmetric matrix S such that Eq. (34) 
holds with Q* orthogonal. Since, S is real and symmetric, it is diagonalizable as 

S = UAU^, 


the eigenvalues in A being real and ordered decreasingly and U an orthogonal matrix. Then, 

Q*UAW + sQ* = e (EhQ* (E°)^ + E^iQ*E°) + a (a^Q* (A°)^ + A^^Q*A°) + F + sQ* 
which can also be written as 

(Q*U)(A + sI)U^ = e (e^Q* (E°)^ + E^,Q*Eo) + a (a^Q* (A^)^ + A^,Q*A°) + sQ* + F. 


(61) 

(62) 

(63) 


In order to conclude the proof, it remains to show that (Q*U)(A + sI)U^ is a singular value decomposition. 
Q*U and U being orthogonal matrices, and A + si being a diagonal matrix, it is sufficient to show that 
A + si has all diagonal positive entries. 

From Eq. 


||A ||2 = ||Q*UAW ||2 = ||e (EhQ* (E°)^ + E^,Q*E°) + « (AhQ* (A°)^ + A^,Q*A°) + F 


(64) 


and 


l|A||2<e 

EhQ* (E°)^ + E^,Q*E° 

T OL 

2 

AhQ* (A0)^ + A^,Q*A° 

<e| 

(||EhI|2||Q1|2 


^ + ||E^,||2||Q*||2||E°||2) 


a 


(iiah||2||q*i|2||(ao) 


|2||Q*||2||A°||2) + ||F||2 


(65) 


<2e||EH||2||E0||2+2a||AH||2||A°||2 + ||F||2 
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by definition of Smin,G in Eq. (30). Denoting by i = 1, • • • the diagonal entries in A, this implies that 


<^111111,0 ^ A^ ^ ^ — f 5 ' ' ' : ^5 (66) 

and, since s > 5min,G: 

A^ H“ s > 0, i = 1, • • • , 

The set of critical points of J'g is included in the set of fixed points of Algorithm 2 and the two sets are 

equal, concluding the proof. 
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